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ABSTRACT 

The solidification of a dilute alloy (bismuth-tin) under Bridgman 
crystal growth conditions is investigated in support of NASA’s 
MEPHISTO space shuttle flight experiment. Computations are 
performed in two-dimensions with a uniform grid. The simulation 
includes the species- concentration, temperature and flow fields, as well 
as conduction in the ampoule. Fully transient simulations have been 
performed; no simplifying steady state approximations are used. Results 
are obtained under microgravity conditions for pure bismuth, and Bi-0.1 
at.%Sn and Bi-1.0 at.%Sn alloys. The concentration dependence of the 
melting temperature is neglected; the solid/liquid interface temperature 
is assumed to be the melting temperature of pure bismuth for all cases 
studied. For the Bi-1 .Oat.%Sn case the results indicate that a secondary 
convective cell, driven by solutal gradients, forms near the interface. The 
magnitude of the velocities in this cell increases with time; this causes 
increasing solute segregation at the solid/liquid interface. 

NOMENCLATURE 

A area 

c p specific heat at constant pressure 

C species concentration, C* / Q* 

D species diffusion coefficient 

E energy 

f volume fraction 

g acceleration due to gravity 

Gr Grashof number, gPiCTH-T^HVv 2 

Gr, solutal Grashof number, gP c Co*H 3 /v 2 

h ampoule thickness (outside radius - inside radius) 

H ampoule diameter; reference length 

i specific enthalpy 

k thermal conductivity 


kp segregation coefficient 

L length of simulation domain 

L a translating zone length 

Le Lewis number, a/D 

Pr Prandtl number, v/a 

q combined convection and diffusion heat flux vector 

q c combined convection and diffusion mass flux vector 

rhsjj coefficient matrix in equation (20) 

RHSy coefficient matrix in equation (15) 
t time 

T temperature 

u,v,w velocities in the x, y, z directions 

Greek Symbols 

a thermal diffusivity 

cty coefficient matrix in equation (20) 

P expansion coefficient 

y solidification front orientation 

A nondimensional temperature, c p * (T- TJ/AH 

AT temperature difference, T H - T c 

AH enthalpy of freezing 

C vorticity 

0 nondimensional temperature, (T-T c )/(Th-T c ) 

Ajj coefficient matrix in equation (20) 

p dynamic viscosity 

p density 

4> nondimensional enthalpy, (i - i,n)/AH 

i {/ stream function 

Subscripts 

0 initial condition 

C cold furnace temperature condition 
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H hot furnace temperature condition 

i, j located at i, j th finite volume center 

I, J located at I, J th finite difference mesh point 

L liquid 

m at solidification front 
mix mixture 

S solid 

w ampoule wall 

Superscripts 
A unit vector 

vector 

* reference or dimensional quantity 

n time step n 

INTRODUCTION 

The synthesis of advanced materials, especially for electronics and 
biomedical applications, demands high-quality crystals. The 
compositional uniformity (and hence the quality) of such crystals can be 
profoundly influenced by the transport phenomena which occur in the 
melt during solidification. The primary transport mechanism causing 
these deleterious effects is natural convection. The low-gravity 
environment of space offers an opportunity to suppress the strength of 
this natural convection. Hence there is a great deal of interest in the 
study of directional solidification of crystals in space. 

The MEPHISTO project (Abbaschian et al. 1992) is a collaborative 
program of space experiments aimed at understanding the fundamental 
processes involved in crystal growth. The space-borne experimental 
apparatus is a Bridgman-type furnace with an isothermal hot zone, an 
isothermal chill zone, and an insulated gradient zone. The MEPHISTO 
furnace contains three ingots of Bi - Sn or Sn- Bi binary alloy inside fused 
silica ampoules with a maximum 6 mm inner diameter and a 1 0 mm outer 
diameter. All three samples are solidified simultaneously, under identical 
thermal conditions. After flight, the samples are extracted and subjected 
to post-flight analysis. Four MEPHISTO space experiments have taken 
place previously; the most recent, MEPHISTO-4, flew in November 
1997. The MEPHISTO-2 and -4 experiments examined the faceted 
solidification of bismuth doped with tin (Bi-Sn). 

The MEPHISTO-2 experiments studied a Bi-0. 1 at.% Sn alloy, with 
the principal aims being: (i) to examine the interfacial morphological 
stability threshold of the faceted solid/liquid interface in the absence of 
convection; (ii) to evaluate the effects of interface kinetics on stability 
during diffusive, fast diffusive and convective transport regimes; (iii) to 
examine Bi growth kinetics as a function of supercooling and (iv) to 
investigate phenomena associated with kinetic roughening (Abbaschian 
et al. 1992). The experimental data from this mission resulted in a 
greater understanding of the dominant role of interface kinetics on 
morphological stability. Stability phenomena were observed that had not 
been previously predicted by theory or measured in terrestrial 
experiments (Abbaschian 1996). 

The MEPHISTO-4 experiment builds on the results of the 
MEPHISTO-2 experiment. MEPHISTO-4 used a Bi- 1.0 at.% Sn alloy. 
Modifications to the experimental apparatus included the addition of a 2 
mm inner diameter capillary tube that extends half-way through the 
process zone. 

The MEPHISTO project includes a program of computational 
modeling of the crystal growth process. In particular, the role of 
convection, which is crucial to a complete understanding of the process, 
is to be investigated. Since accurate experimental determination of 


convection in metallic melts is very difficult to achieve, due to the 
opacity and chemical reactivity of the melts, convective levels are 
determined numerically. Furthermore, the computational models 
themselves are to be improved by a process involving prediction of, and 
comparison with, the experimental results. The aim of this procedure is 
to develop effective fully three-dimensional computer simulations of 
fluid flow related effects (Abbaschian 1996). Previously, three- 
dimensional solutions for Bridgman growth have been limited to steady- 
state growth in succinonitrile, a widely used transparent phase change 
material with properties analogous to metallic materials (Yao and de 
Groh 1993, de Groh and Yao 1994). 

As part of pre-flight computational analysis in support of 
MEPHISTO-4, convection effects at microgravity levels and the 
influence of the 2 mm diameter growth capillary have been modeled by 
means of a transient, 2D FIDAP finite-element model (Yao et al. 1997). 
A fixed-grid approach was used, with the enthalpy method being 
employed to model the phase change. Temporal averaging was used for 
the apparent heat capacity in the discretized equations. Due to 
computational difficulties introduced by the small partition coefficient for 
Bi-Sn, the presence of solute was ignored. 

Preliminary scaling arguments of the convection levels by de Groh 
and Nelson (1994) imply that solutal convection effects on solute 
segregation may be significant. However, thus far, it has not been 
possible to include solutal convection, or even a passive solute, into 
numerical simulations involving phase change for MEPHISTO-4. This 
is due to the difficulties concerning convergence with front-tracking 
methods as well as those imposed by the low value of the partition 
coefficient for Bi-Sn (Yao et al. 1995, Yao et al. 1997, Simpson et al. 
1998). 

Many simulations of Bridgman crystal growth processes, both under 
terrestrial and microgravity conditions, are available in the literature. 
The majority of these simulations can be classified as pseudo steady state 
models (PSSM). The key assumption in such models is that a “steady 
state” mode of alloy solidification exists, i.e. the concentration of the 
dopant in the solid which forms at the interface is equal to the initial 
dopant concentration in the liquid (Kurz and Fisher 1989). Such models 
vary in complexity from simple 2D analyses that consider the interface 
to be flat (Alexander et al. 1989, Yao et al. 1995, Simpson et al. 1998) 
to much more complex formulations that are able to handle interface 
curvature and wall conduction (Adomato and Brown 1987) and fully 
three dimensional simulations (Yao and de Groh 1993, Liang and Lan 
1996). However, because of the low partition coefficient for Bi-Sn 
alloys, a “steady state” mode of solidification is never achieved during 
the MEPHISTO experiments. Thus pseudo steady state models are not 
appropriate; fully transient simulations need to be developed in order to 
faithfully model this process. 

The computational modeling presented in this paper is intended to 
examine the effects of thermosolutal natural convection on the 
MEPHISTO-4 solidification experiments. This will be achieved using 
a fully transient 2D model, which includes most of the effects of binary 
alloy solidification and convection driven by both thermal and solutal 
gradients. Conduction through the ampoule wall will be considered, 
while the capillary tube and the dependence of melting temperature on 
concentration will be ignored. For the purposes of our model, we assume 
that the liquid/solid interface remains distinct like that for a pure material. 

MATHEMATICAL FORMULATION 

The problem under consideration is the directional solidification of 
a binary alloy by the Bridgman process, as shown schematically in Fig. 1 . 
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Fig. 1. Schematic of the Bridgman crystal growth process and furnace 
temperature profile. 


The gravity vector is perpendicular to the furnace axis (horizontal 
Bridgman growth configuration). The melt region is considered to be a 
viscous Newtonian fluid subject to thermosolutal convection. 
Thermophysical properties are considered as constant but distinct for the 
solid and liquid phases. Density variations are considered to be subject 
to the Boussinesq approximation. The governing equations for the 
velocity field are the vorticity and vector potential equations which, in 
nondimensionalized form, are (Mallinson and de Vahl Davis 1973): 

+Vx(£xO) = -Gr(Vx0g) -Gr s (VxCg) + PrV 2 £ 0) 
dt 

V 2 ijr = -C. W 

In these equations, the definition of vorticity and the relationship 
between velocity and vector potential are 

( = -Vxu, u = Vxijf (?) 

In the nondimensionalization, the ampoule inside diameter H, is selected 
to be the reference length. The characteristic time and velocity become 
f (= H 2 /a) and v* (= H/f = aJ H). 

The solution of equation (1) subject to equation (2) constitutes the 
vorticity- vector potential formulation for natural convection. For a non- 
slip wall, the boundary conditions to be imposed on vorticity and vector 
potential (Mallinson and de Vahl Davis 1973) are that the velocity is 
zero, as are the tangential derivatives of its components. This leads to the 
following conditions for vorticity 


x =0 
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(4) 


The vector potential at a plane, impermeable surface is normal to the 
surface and its gradient is zero: 



x = 0 i|f 2 = 0 

i |/ 3 = 0 


( 5 ) 


The boundary conditions for the other surfaces result from similar 
considerations (Mallinson and de Vahl Davis 1973). As an initial 
condition, the fluid must be quiescent 

t = 0 <p = C = 0. ( 6 > 

The nondimensionalized form of the governing equation for the 
conservation of energy is 

— +V-(u0)=V 2 0 (7) 

at 


It will be seen below that an integral form of this equation is suitable for 
calculating phase-change. An initial temperature equal to the hot wall 
temperature is applied throughout the flow field. The boundary 
conditions at the x = 0 and x = L walls are shown in Fig. 1 . 

The thermal boundary conditions along the y = -h and y = H+h walls 
required for the Bridgman process need closer examination, since they 
are a function of time. A schematic of the temperature imposed at these 
boundaries, which approximates the furnace, is shown in Fig. I. There 
is a translating zone (considered an “adiabatic” zone if the temperature 
profile is unknown) between the hot and cold regions of the furnace, in 
which the temperature linearly increases from the cold furnace 
temperature to the hot furnace temperature. The melting temperature of 
the material occurs somewhere within this zone, which translates with 
time at a constant x-velocity, known as the translation velocity , u l . This 
is what facilitates the directional growth of the crystal. Defining the x 
location where the translating zone meets the cold furnace temperature 
zone to be at x A (t), the boundary condition for temperature may be 
expressed as 

T c , for x < x A (t) 


at y = -h, H+h: T 


x -x A (t) 
L a 


, for x A (t) as [x A (t)+L A ] (8) 


T h , for x A (t) < x 

In principle, the solution of the energy equation (7) coupled with the 
solution of the vorticity- vector potential equations (1), (2) and (3) would 
yield the temperature and velocity distribution throughout the simulation 
domain. However, the problem of modeling the physics of the 
propagation of the solidification front and determining its location 
remains to be addressed. We choose to do this by the enthalpy method 
(Alexiades and Solomon 1993). The constitutive equation required for 
this method is the integral equation describing the conservation of energy 
in an arbitrary control volume: 


i*At 


- q • n dA dt 


(9) 


Unlike the differential form of the energy equation (7), equation (9) is 
valid for any control volume in the solution domain, including control 
volumes which are intersected by the interface. 

The nondimensionalized equation for the conservation of solute 
throughout the computational domain is 


— + V-(uC) = — V 2 C (10) 

3t Le 


This equation is analogous to the energy equation. We impose an initial 
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solute concentration throughout the solution domain. At the boundaries 
we demand that there be no flux of solute: no solute may exit the solution 
domain. Thus, 

1 = 0 C = C 0 

x = 0, L dCJdx = 0 (11) 

y = 0, H dC/dy = 0 

Again, in principle, solution of equation (10) along with energy 
equation (7) and vorticity- vector potential equations (1), (2) and (3) all 
subject to the relevant boundary and initial conditions are enough to 
determine the solute, temperature and velocity values throughout the 
solution domain. However, the more general problem involving phase 
change demands that the thermodynamics of solute redistribution be 
addressed. The methodology for this is described in the next section. 
The constitutive equation for this method uses the integral form of 
equation (10), 


l*Al l*Ai 

(s(/ c -* dv ) d, = // 


- q c • n dA dt (12) 


NUMERICAL ANALYSIS 

The computational domain is primarily discretized using regularly 
spaced finite difference mesh points. Superimposed on this grid are finite 
volumes which are required for the enthalpy method. The finite- volume 
centers are staggered with respect to the finite-difference mesh point 
locations. Vorticity, velocity and vector potential are calculated at the 
finite difference mesh points. Temperature, solute concentration and 
enthalpy are evaluated at the control volume centers as a result of the 
solution of equations (9) and (12). 

The discretization and the accompanying solution scheme are 
explained in Simpson and Garimella (1997). The approach used here is 
modified from a program (Timchenko et al. 1997) written for the solution 
of natural convection in a rectangular cavity. The essential details are 
that the vorticity transport equation (1) and the vector potential equation 
(2) are discretized in space using the regular finite difference mesh. The 
discretized equation for the vorticity transport equation is solved using 
an Alternating Direction Implicit (ADI) scheme (Samarskii and Andreyev 
1963). The discretized equation for vector potential is solved using the 
conjugate gradient method. This is a semi-iterative method used to solve 
specific systems of linear equations. The method is not discussed here, 
except to note that the values for vector potential at the previous time 
step are used as the initial guess values. Once the values of vector 
potential are known, the nodal velocities can be determined from 
equation (3). 

It should be noted that while governing equations (1) and (2) are 
presented in three dimensions and the code is capable of solving for the 
flow field in 3D, the simulations in this paper are restricted to two 
dimensions. This is accomplished by solving only the z-component of 
vorticity and vector potential and thus yielding the requisite u- and v- 
velocity components. 

The discretization for the boundary conditions on vorticity (4) and 
vector potential (5) at a wall aligned with the Cartesian grid are easy to 
realize. However, the boundary conditions at the arbitrarily oriented 
interface require special handling. The slope of the interface is found 
using a Hirt and Nichols (1981) type front reconstruction; the boundary 
conditions may then be applied (following Raw and Lee 1991, Roache 
1976) as: 

4*3 =0 C 3 = (cos' 2 Y)5 2 vj/ 3 /0y 2 = (sin' 2 Y)5 2 ^ 3 /^x 2 


The slope of the interface is used to determine which of the two 
formulations for vorticity at the boundary is used. If the line is more 
horizontal than vertical then the first formulation is used and the converse 
applies. 

The integral energy equation (9) is discretized using the finite 
volume mesh. An upwind scheme is incorporated for the treatment of 
convective heat fluxes. This discretization may then be written as 


— -AxAy = 
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in which [ A, B ] is the maximum of A and B. The velocities at the finite 
volume faces are the interpolated values from the finite-difference mesh 
points. 

The modified enthalpy method is the solution scheme for the energy 
equation. The thrust of the method is to use equation (9) to determine the 
cell temperature and enthalpy. The complete methodology for the 
solution is explained in Simpson and Garimella (1997); an explanation 
of the modifications necessary to account for the presence of the ampoule 
wall is provided in Garimella and Simpson (1998). Essentially, equation 
(14) is advanced forward in time using Gauss-Seidel iteration with 
successive over relaxation. For this method, equation (14) reduces to 

d)?* 1 +C ? A?* 1 = RHS ? O 5 ) 

t g g g m 

in which p denotes the inner iteration number and the RHS matrix is fully 
explicit at iteration p. The values at time step n are used as the first 
approximation. In the melt region or in the solid, the relationship 
between temperature and enthalpy gives us 


(J ) p : 1 = 


RHS L P 

1 ±-, if RHS^ <; 0 


C • + c 


0, 

RHsg+cg 
<^c pI ’ 


if 0 < RHSjj < 1 


(16) 


if RHS^ ;> 


with corresponding results for the temperature (X). The scheme is 
modified in the manner explained in Garimella and Simpson (1998) so 
that it is also valid for the ampoule wall region and furnishes the 
temperatures in that region without any special handling. Iterations 
proceed until convergence, at which time the new values at time step n+ 1 
for temperature and enthalpy are declared to be those found at 
convergence. Convergence is assessed by examining the change in the 
solution and the change in the solution residual with a tolerance of 10' 6 
(Garimella and Simpson 1998, Barrett et al. 1994). In cells which 
contain the interface, the cell enthalpies are equivalent to the cell liquid 
fractions by virtue of our nondimensionalization scheme. This allows the 
front orientation to be determined using Hirt and Nichols (1981). The 
front orientations are required for the calculation of vorticity at the 
interface. 

The discretized analog of equation (12) is 
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The assumptions made in arriving at this equation are that the densities 
of the liquid and solid phases are constant and equal, so that a simple 
mixture rule applies, and that there is no diffusion in the solid (except in 
the cells that contain the interface, in which the equilibrium lever rule is 
assumed); the concentration at which the solute first solidifies is the 
concentration at which that portion of solid remains for all time. 

Solution of this finite volume equation is explained in Simpson and 
Garimella (1998). For the purposes of our model, we assume that there 
is no effect of concentration on the melting temperature and the 
liquid/solid interface remains distinct like that for a pure material. The 
significant result of this assumption is that the enthalpy method is 
decoupled from the concentration equation and is the only mechanism 
necessary to determine the location of the solidification front. 

The LHS of equation (17) may be expressed in terms of liquid solute 
concentration, C L as 




At 


At 


AxAy 


(18) 


in which 

<j = [< k p" , > f " ,+(2 ‘ k P> f ».*) 

and this is known at time step n since the liquid fractions at n+1 are 
recovered directly from the solution of the energy equation. The solute 
concentration in the solid portion of a node may be found from 
employing the relation 


C 


n ♦ 1 
S 


Cs°fs"^ p C L n (f L n>1 -f L n ) 


c-n ♦ 1 


(19) 


The discretized equation (17), subject to (18) is solved by Gauss- 
Seidel iteration in a manner analogous to the energy equation, i.e., 

k^Jc^'-rhs^ (20) 

The solution scheme for simplified, fully transient solidification in 
two-dimensional Bridgman crystal growth subject to thermosolutal 
convection has now been fully specified. 


RESULTS AND DISCUSSION 

Pure Bismuth 

Simulations for the Bridgman crystal growth of pure bismuth were 
performed. Thermophysical properties from Yao et al. (1995) were used. 


The cold and hot furnace temperatures were T c = 50°C and T„ = 700°C 
respectively. Liquid properties were used as the reference properties for 
the nondimensionalization scheme. The properties for solid and liquid 
bismuth were evaluated at the mean solid and the mean liquid 
temperatures of 160.7°C and 485.7 °C, respectively. The reference 
density was considered to be equal for both phases. The property values 
were: k, = 6.872 W/mK, k, = 14.66 W/mK, ft = 10,070 kg/m 3 , = 132.6 
J/kgK, Cpj = 135.3 J/kgK, AH = 52.3 kJ/kg, p = 1.240 x 10 3 Ns/m and T m 
= 271.3°C. The thermal expansivity of the liquid was taken to be p x = 
- 1 .25 x 10 4 K 1 . The properties used for the alumina ampoule wall were 
k w = 2.10 W/mK, p w = 2,020 kg/m 3 , c pw = 92.10 J/kgK (In the 
MEPHISTO flight experiments, a fused silica ampoule is used, which has 
a higher value for specific heat. This is not expected to have a significant 
impact on the results. In future work, fused silica properties will be 
employed). The dimensions used to define the domain were H = 6 mm 
and L = 73 mm. The ampoule wall thickness, h, was 2 mm. The length 
of the insulated translating zone was L A = 25 mm, with a translation 
velocity of = 3.38 pm/s. The gravity level was 10 pg, with gravity 
acting in the negative y direction. The important nondimensional 
quantities from these values are Gr = 1 13.6 and Pr = 0.01 144. In order 
to start the crystal growth simulations, the following procedure was 
carried out in each simulation. The initial position of the translating zone 
was flush with the x = 0 wall. This zone is immobilized for the first 

5.000 time steps. During this time, the velocity and concentration field 
solution schemes are switched off while solid rapidly chills in the portion 
of the translating zone which is lower than the melting temperature. 
After this, simulations proceed with the entire solution scheme enabled, 
and the insulated zone moving at the translation velocity u,. 

The computational domain was discretized by means of a regularly 
spaced 200 x 28 mesh, with 20 of the 28 y-direction cells being in the 
sample and the remaining 8 in the ampoule walls. This mesh was 
selected as a result of numerical experiments which indicated its 
adequacy in resolving solute gradients in the liquid near the interface. 
The temporal step size was At = 6.25 x 10' 2 . The program was written 
in ANSI standard FORTRAN 77 and implemented on a DEC 3000/700 
(225 MHZ, SPECfp95 5.71). The total CPU requirement for 61,000 time 
steps was 4.33 hours. Note the efficiency of the method even with the 
fine mesh. 

Figure 2 (a) shows a plot of velocity vectors and isotherms after 

19.000 time steps (0.81 hours) have elapsed. Isotherms are shown in the 
ampoule wall as well as in the liquid and solid bismuth. The thick line 
at x = 3.15 represents the solid/liquid interface. The dominant feature of 
this plot is the counter-clockwise convective cell in the translating zone 
(3.15 < x < 7.8). For clarity, velocity vectors are plotted on every other 
mesh point in the x-direction in this figure. The maximum 
(nondimensional) velocity in this convective cell is 1.936 x 10 3 at 
(4.625, 0.2). The ratio of this velocity to the translation velocity, u t , is 
1.08, which compares well with the value of 1.12 found using a finite- 
element, variable-property simulation of the same process by Yao et al. 
(1997). Velocities in the negative y-direction at the interface are 
constrained and concentrated by the presence of the solid whereas 
velocities in the positive y-direction at the opposite (hot) end of the 
translating zone are much more diffuse. 

The other main feature of this plot is the isotherms throughout the 
solution domain. On the outside edge of the ampoule, the (imposed) 
temperature profile is a linear ramp function within the translating zone. 
This is easily witnessed by the regularly spaced isotherms on the outer 
edge. In the ampoule region, the isotherms are dramatically distorted. 
Since the ampoule material has a low thermal conductivity, the thermal 
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field on the inside of the ampoule wall differs from that imposed on the 
outside. For xs 4.4, the temperature in the ampoule is greater than the 
applied temperature (note the isotherms distort sharply to the left). For 
4.4, the converse is true. In the bismuth region the isotherms exhibit 
the same trend, and have a gentle, crescent-shaped curvature. The 
isotherms and the interface appear to be symmetric about the centerline 
(y = 0.5) and thus have not been discemibly influenced by convective 
transport in the melt. 

1 2 9 4 5 6 7 5 

■ 0.4554 - 0.2555 - 0.0456 0.1643 0-3743 0.5542 0 . 7*41 1.0040 



Fig. 2. Velocity vectors and isotherms for the Bridgman growth of pure 
bismuth at (a) 19,000 time steps (0.81 hours), (b) 37,000 time steps (1.63 
hours) and (c) 55,000 time steps (2.90 hours). The solid line indicates 
the location of the solid/liquid interface. The velocity vectors, shown at 
every other location in the x-direction, indicate a single, counter- 
clockwise convective cell. 

Figures 2 (b) and 2 (c) are plots of velocity vectors and isotherms at 
37,000 (1.63 hours) and 55,000 time steps (2.90 hours), respectively. 
The progression of the translating zone (and thus of the solidification 
front) is obvious from these figures. After 1.63 hours, the front is at x « 
5.2, and the nature of the convective motions is unchanged. The 
magnitude of the maximum velocity is identical to that at the previous 
time (Fig. 2a), and its new location is (6.6875, 0.2). At a still later time 
(2.90 hours), when the front is at x « 7.4, the velocity increases slightly 
to 1 .943 x 10" 3 and is at (8.8125, 0.2). This indicates that end effects are 
coming into play, and so, results after this point are not considered. 

Figure 3 is a trace of temperatures at three different y-locations, 1 .63 
hours into the process (Fig. 2b). The solid line (y = 1.2917) is the 
applied thermal boundary condition. The y = 1.0 trace is along the inside 
of the ampoule wall. These two traces are significantly different, 
indicating the need to include the ampoule wall in the simulation. The 
third trace (y = 0.5) is the temperature along the centerline of the domain. 
Note the change in slope at the solidification front (at x = 5.2). This is 
due to the change in material properties in the liquid and solid phases. 
These results agree well with those found in Yao et al. (1997), both in 
terms of the temperature values and the slope change due to property 
variation at the interface. 



Fig. 3. Temperature traces at y = 1.2917 (the outermost cell in the 
ampoule), y = 1.0 (along the inside of the ampoule wall) and y = 0.5 (the 
centerline), for the pure bismuth case after 1 .63 hours. The results show 
the dramatic influence of the ampoule on the thermal field applied to the 
bismuth sample. 

Bi - 0.1 aL% Sn alloy (MEPHISTO-2) 

For simulating the Bridgman crystal growth of a Bi- 0.1 at.% Sn 
alloy, the thermophysical properties used were the same as for the 
simulation of pure bismuth described above. In addition, the value for 
diffusivity of liquid Sn in liquid Bi was D = 6.22 x 10 9 m 2 /s, and the 
partition coefficient was kp = 0.029. The value of solutal expansion 
coefficient was taken to be p c = -0.305 (at.%)' 1 . The diffusivity value was 
calculated at the mean melt temperature. The liquid composition was 
taken to be at a uniform value of 0. 1 at.% as the initial condition. These 
values result in dimensionless parameters Gr, = 0.4262 and Le = 1730. 
Simulations were performed using the same computer, spatial 
discretization and time step as for the pure bismuth case. The CPU 
burden for 60,000 time steps was 4.98 hours. 

Velocity vectors and isotherms after 0.81 hours are shown for this 
alloy in Fig. 4 (a). The nature of the convective flow is similar to that for 
pure bismuth (Fig. 2a). A single, counter-clockwise, thermally driven 
convective cell dominates the velocity plot. The maximum u-velocity is 
slightly less than for pure bismuth, and is at the same location. However, 
the magnitude of the maximum v-velocity is 9.874 x 10 4 at (3.375, 0.5) 
acting in the negative y direction, which is lower in magnitude than for 
pure bismuth (1.047 x 10 3 ). Solute rejected at the interface acts to 
oppose the thermally driven convective motion, resulting in retarded 
velocities near the interface. The thermal field is not distorted by the 
action of convection. Figure 4 (b) is a plot of velocity vectors and 
isotherms after 1.63 hours. Again the convective motion and temperature 
field is similar to the pure bismuth case (Fig. 2b). The magnitude of the 
maximum v-velocity has decreased to 9.527 x 10' 4 acting downward at 
(5.5, 0.5) as a result of continuing solute rejection at the interface. At a 
later time (2.90 hours. Fig. 4c), the maximum v-velocity is 9.353 x 10 4 
(at 7.6875, 0.5) confirming the trend of solute rejection opposing the 
convective motion due to the applied temperature gradient. 
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Fig. 4. Velocity vectors and isotherms for the Bridgman growth of 
Bi-0.1 at.% Sn at (a) 0.81 hours, (b) 1.63 hours, and (c) 2.90 hours. The 
solid line indicates the location of the solid/liquid interface. A single, 
counter-clockwise convective cell is present in all cases. 


Solute concentrations in the liquid at the same three times as in Fig. 
4 are illustrated in Fig. 5. The solute build-up near the interface, and the 
exponential profile characteristic of binary-alloy solute rejection, are 
clearly evident. The increased solute concentration levels at small y 
result from the convective flow sweeping solute “down” the interface in 
the direction of decreasing y and away from the interface in the direction 
of increasing x (refer to Fig. 4). Closer to the interface, the solute 
concentration exhibits a minimum value near y = 0.7. The minimum 
value is found at this location for two reasons. The curvature of the 
interface causes the value at the centerline to be lower than that found at 
the outside edges. In the absence of convection the minimum value 
would be expected to occur at y = 0.5. In addition, the incoming hot, 
solute-poor liquid in the upper portion of the domain causes this 
minimum value to be offset to y = 0.7. Note also that the flow near the 
top is retarded due to the presence of the non-slip wall; this results in 
increased solute concentration levels there. 


y 


y 


y 



U U iJO U 10 



M U 7JB 



7 S 10 Li 9J0 


X 


(a) 


(b) 


(c) 


Fig. 5. Liquid solute concentration values for Bi-0.1 at.% Sn at (a) 0.81 
hours, (b) 1.63 hours, and (c) 2.90 hours. The profiles are a result of 
combined diffusive and convective processes. The maximum 
concentration increases for each case, indicating solute build-up. 

In Fig. 5 (a), the extreme values for solute concentration at the first 
mesh points adjacent to the interface (x = 3. 1 563) are 0.4207 at.% at y = 
0.175 and 0.3718 at.% at y = 0.675. The front is approximately at x = 
5.2 after 1.63 hours have elapsed (Fig. 5b). Solute buildup is seen to 
continue as the front progresses, with the extreme values of liquid solute 
concentration near the interface (x = 5.2813) being 0.5894 at.% at y = 
0.875 and 0.5347 at.% at y = 0.625 respectively. Although the maximum 
value is now near the top of the domain (whereas it was near the bottom 
in Fig. 5a), the value for solute concentration near the bottom (y = 0.275) 
is still as high as 0.5887 at.%, indicating that solute segregation in the y- 
direction is not monotonic. In Fig. 5 (c) the front is at approximately x 
= 7.4, and the extreme solute concentrations near the interface (x = 
7.4063) are 0.7105 at.% at y = 0.125 and 0.6429 at.% at y = 0.625. Note 
that these concentration values will not cause a significant change in the 
melting temperature of the alloy as predicted by the Bi-Sn binary phase 
diagram. Thus the assumption of a constant melting temperature is valid. 
Figure 5 shows that the solute levels near the interface increase with time. 
It appears that this solute build-up is responsible for a corresponding 
decrease in the convective velocities near the interface with time (refer 
to the preceding discussion of the velocity field). The presence of solute 
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near the interface acts to make the fluid less dense. This opposes the 
action of the thermal field which increases the fluid density near the 
interface. The net effect is that the thermally driven convection near the 
interface is retarded by this opposing potential. 

For the low-concentration alloy, the liquid solute concentration 
profiles are governed by solute rejection at the interface and by thermal 
convection only. The action of solutal convection is limited at these low 
concentrations. The velocity field indicates that thermal convection 
decreases slightly with time. 

Bi -1.0 at.% Sn alloy (MEPHISTO-4) 

A richer alloy was simulated next. The liquid composition was 
taken to be at a uniform value of 1.0 at.% Sn as the initial condition 
(yielding Gr, = 4.262). Unlike the previous simulations, computations 
were performed using a Silicon Graphics Indigo 2 workstation (250 MHZ 
CPU, SPECfp92 177.5). The CPU burden for 61,000 time steps was 
8.27 hours. 
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Velocity vectors and isotherms after 0.81 hours for this alloy are 
shown in Fig. 6 (a). The velocity vectors indicate a primary convective 
cell rotating in a counter-clockwise manner in the translating zone. 
However, the velocities differ from those shown in Figs. 2 (a) and 4 (a); 
the convective velocities in the region near the interface are significantly 
smaller in Fig. 6 (a). In fact, there are some positive v- velocities in the 
region near the interface at y =* 0.9. Since these positive velocities are 
only present at a few mesh points, an organized convective cell driven by 
solute gradients is not discerned; however, a v-velocity contour plot (not 
shown) revealed the presence of a small convective cell. The maximum 
u-velocity for this case is +1.885 x 10' 3 at (4.8125, 0.2), which is lower 
than the value of 1.936 x 10‘ 3 found for pure bismuth. The maximum v- 
velocity magnitude is 7.702 x 10‘ 4 acting downward at (3.6875, 0.5), 
which is also significantly lower than the value of 1.047 x 10‘ 3 for pure 
bismuth, or the 9.874 x 10' 4 value found for Bi-0.1 at.% Sn. Higher 
levels of solute rejection due to the higher initial concentration of the 
melt are responsible for these effects. 
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Fig. 6. Velocity vectors and isotherms for the Bridgman growth of 
Bi- 1.0 at.% Sn at (a) 0.81 hours, (b) 1.63 hours, and (c) 2.90 hours. The 
solid line indicates the location of the solid/liquid interface. A primary 
counter-clockwise convective cell, driven by thermal gradients, is evident 
in all three panels. A secondary, clockwise cell, driven by solutal 
gradients, develops with time. 


Fig. 7. Liquid solute concentration values for Bi- 1 .0 at.% Sn at (a) 0.8 1 
hours, (b) 1.63 hours, and (c) 2.90 hours. The profiles are a result of 
combined diffusive and convective processes. The maximum 
concentration increases for each case, as does the segregation in the y- 
direction. 
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At a later time (1.63 hours), the front is at x = 5.4, as shown in Fig. 
6 (b). A secondary, clockwise convective cell driven by solute gradients 
has formed adjacent to the interface. The two-cell convective motion is 
in contrast to that observed in Figs. 2 (b) and 4 (b). The v- velocities near 
the interface are positive along the entire interface height. The maximum 
v-velocity in the secondary convective cell near the interface at (5.375, 
0.7) is +2.187 x 10'\ and the maximum u-velocity in the domain is 
+1.858 x 10 3 at (7.0, 0.2). Later in the growth process (2.90 hours), the 
front has advanced to x = 7.4, and the secondary convective cell has 
increased in size and strength as solute - the driving force for this cell - 
continues to build up at the interface. The maximum u-velocity in the 
domain is +1.852 x 10' 3 at (9.1875, 0.2), while the maximum v-velocity 
in the secondary cell (and in the domain) is +4.441 x 10' 4 at (7.5, 0.6). 

Plots of the solute concentration corresponding to Fig. 6 are shown 
in Fig. 7. As for the Bi-0.1 at.% Sn case, this plot for the richer alloy 
exhibits solute build-up near the interface, and the exponential profile 
characteristic of binary-alloy solute rejection. In Fig. 7 (a), close to the 
interface, the level of solute concentration increases with increasing y. 
Further from the solute interface, but still within the solute boundary 
layer, concentration increases with decreasing y, for a given x value. 
This is due to the localized convection in the positive y direction near the 
interface and the bulk thermally driven convection acting elsewhere (refer 
to Fig. 6a). At x = 3.1562, the maximum concentration is 4.559 at.% at 
y = 0.875, while the minimum value is 3.663 at.% at y = 0.175. Interface 
curvature effects are also evident, with the profile being preferentially 
rich away from the centerline (y = 0.5). 

As time progresses (Fig. 7b), the liquid becomes more solute-rich, 
with the maximum and minimum values near the interface (x = 5.2812) 
being 7.205 at.% at y = 0.825 and 4.747 at.% at y = 0.175, respectively. 
The additional solute rejected and stronger solutal convection (see Fig. 
6b) have led to higher levels of solute segregation in the y-direction when 
compared to Fig. 7 (a). The higher levels of solute rejection and solutal 
convection in the liquid will cause higher levels of segregation in the 
solid. Figure 7 (c) is a plot of liquid solute concentration after 2.90 
hours. The general character of the profile has changed slightly. The 
existence of the secondary cell has decreased the level of convection in 
the negative y-direction in the outer regions of the boundary layer. Close 
to the interface, the strong positive vertical velocities have led to 
increased solute build-up at large y. The maximum and minimum values 
of interfacial liquid solute concentration (x = 7.4062) are 9.682 at.% at 
y = 0.875 and 5.423 at.% at y = 0.175, respectively. Note that 
concentration values of this magnitude would be sufficient to lower the 
melting temperature of the alloy by a significant amount (see Simpson et 
al. 1998). This implies that the assumption that the melting temperature 
remains constant is not valid. The likely effect of concentration- 
dependent melting temperature would be to cause the interface to be 
thicker near the bottom of the domain. This effect may alter the 
concentration and velocity fields and will be addressed in future work. 

For the higher alloy concentrations, solutal convection plays a much 
larger role, as can be seen in the concentration profiles discussed above. 
As the growth proceeds, the level of solutal convection grows, with a 
corresponding change in the solute distribution. The maximum 
concentration increases much more rapidly than the minimum value. 

CONCLUSIONS 

A series of fully transient simulations of horizontal Bridgman crystal 
growth in microgravity conditions have been performed. The pure 
bismuth simulation was found to agree well with a similar simulation of 
the process (Yao et al. 1997), in terms of the convection level and the 


presence of slight interface curvature. 

For the dilute alloy simulation (Bi-0.1 at.% Sn), a single dominant 
counter-clockwise rotating convective cell is present for the entire 
duration of the process. As time proceeds, solute is rejected at the 
interface and the level of solute near the interface increases. Due to the 
presence of convection, radial solute segregation occurs with 
preferentially higher values at the bottom of the domain. It is unclear if 
the level of segregation increases with time. The maximum values of 
concentration are small and so the assumption of constant alloy melting 
temperature is realistic. 

For the richer alloy (Bi-1.0 at.% Sn), the convective field is much 
more complex. Higher levels of solute rejection at the interface cause 
higher levels of solutal convection. Initially, a single, thermally driven, 
counter-clockwise rotating cell is present. As time proceeds and solute 
accumulates at the interface, a secondary, solute driven clockwise 
rotating cell develops. The effect of this convective pattern is to yield 
significant levels of segregation in the solute near the interface (and 
hence in the solid). This segregation is such that the values near the top 
of the domain are maximum. The higher values of concentration near the 
interface would result in a significant change in the melting temperature 
for the alloy. Thus the assumption that the melting temperature is 
constant is not valid for this case. Investigating the effects of 
concentration-dependent melting temperature will be addressed in future 
work. A comparison of results with existing microgravity experimental 
data and a further series of simulations for Bridgman growth of 
succinonitrile (Yao and de Groh, 1993) under terrestrial and microgravity 
conditions are also planned as future work. 
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